home *** CD-ROM | disk | FTP | other *** search
/ Mac Mania 5 / MacMania 5.toast / / Tools&Utilities / Plotfoil 3.2 / naca-1.0 / naca4.f < prev    next >
Text File  |  1995-09-13  |  9KB  |  377 lines
  1. c
  2. c-----------------------------------------------------------------------------
  3. c
  4. c  Naca4.f -- contains routines to generate the NACA four and five
  5. c             digit series foils, both in their original form, and
  6. c             with the modified thickness profiles.
  7. c
  8. c  Written By: S.E.Norris
  9. c
  10. c  norris@cfd.mech.unsw.edu.au
  11. c
  12. c  $RCSfile: naca4.f,v $
  13. c  $Author: norris $
  14. c  $Revision: 1.3 $
  15. c  $Date: 1995/08/31 11:05:52 $
  16. c
  17. c  $Log: naca4.f,v $
  18. c  Revision 1.3  1995/08/31  11:05:52  norris
  19. c  *** empty log message ***
  20. c
  21. c  Revision 1.2  1995/08/31  06:11:57  norris
  22. c  Fixed divide error that was stuffing up on the Cray.
  23. c
  24. c  Revision 1.1  1995/08/30  09:59:47  norris
  25. c  Initial revision
  26. c
  27. c-----------------------------------------------------------------------------
  28. c
  29.       SUBROUTINE Naca4( x,npl,npmx,naca,inaca,scle,modified )
  30. c
  31.       IMPLICIT none
  32.       INTEGER npl,npmx,inaca
  33.       REAL x(3,npmx),scle
  34.       CHARACTER naca*(*)
  35.       LOGICAL modified
  36. c
  37. c  This routine calculates the ordinates of a NACA four digit foil.
  38. c
  39.       INTEGER i,j
  40.       INTEGER npt,npf
  41.       REAL ndata(5),cf(3),ca(3),ct(5),dcf(2),dca(2)
  42.       REAL a(5),d(4)
  43.       REAL dx,dy,dydx,t,y
  44.       REAL Poly
  45.       EXTERNAL Poly
  46. c
  47.       DATA  ct / 0.0 ,-0.126 ,-0.3537 , 0.2843 ,-0.1015 /
  48. c
  49. c
  50.       npt  = npl+1
  51.       npf  = npl/2
  52. c
  53. c     Set the x values according to a cosine distribution.
  54. c
  55.       call cosdn( x,npl,npmx )
  56. c
  57. c     Calculate the parameters describing the foil, from the given NACA
  58. c     number. Then calculate the polynomial coeffs for the camber.
  59. c
  60.       call Naca4b( naca,inaca,ndata,modified )
  61.       call Naca4c( cf,dcf,ndata(1),ndata(2) )
  62.       call Naca4c( ca,dca,ndata(1),(1.0 - ndata(2)) )
  63.       if (modified) then
  64.      call Nacam( a,d,ndata(4) )
  65.       endif
  66. c
  67. c     Calculate 'y' values for each value of 'x'.
  68. c
  69.       x(2,1)     = 0.0
  70.       x(2,npt)   = 0.0
  71.       x(2,npf+1) = 0.0
  72.       do i = 2,npf
  73.      j = npt+1-i
  74. c
  75. c        Calculate the camber line.
  76. c
  77.      if (x(1,i).le.ndata(2)) then
  78.         y    = Poly( x(1,i),2,cf  )
  79.         dydx = Poly( x(1,i),1,dcf )
  80.      else
  81.         y    = Poly(( 1.0 - x(1,i) ),2,ca  )
  82.         dydx = Poly(( 1.0 - x(1,i) ),1,dca )
  83.      endif
  84. c
  85. c        Calculate thickness.
  86. c
  87.      if (modified) then
  88.         if (x(1,i).le.ndata(5)) then
  89.            a(1) = a(5)*sqrt( x(1,i) )
  90.            t = ndata(3)*Poly( x(1,i),3,a )
  91.         else
  92.            t = ndata(3)*Poly( (1.0-x(1,i)),3,d )
  93.         endif
  94.      else
  95.         ct(1) = 0.2969*sqrt( x(1,i) )
  96.         t = 5.0*ndata(3)*Poly( x(1,i),4,ct )
  97.      endif
  98. c
  99. c        Finally, calculate node positions upon the foil surface.
  100. c
  101.      dy = sqrt(( t**2 )/( 1.0+( dydx**2 )))
  102.      dx = dy * dydx
  103. c
  104.      x(1,i) = x(1,i) + dx
  105.      x(1,j) = x(1,i) - 2.0*dx
  106.      x(2,i) = y - dy
  107.      x(2,j) = y + dy
  108.       enddo
  109. c
  110. c     Rescale dimensions if nessisary.
  111. c
  112.       call FoilScale( x,npt,scle )
  113. c
  114.       return
  115.       END
  116. c
  117. c-----------------------------------------------------------------------------
  118. c
  119.       SUBROUTINE Naca4b( naca,inaca,ndata,modified )
  120. c
  121.       IMPLICIT none
  122.       INTEGER inaca
  123.       REAL ndata(5)
  124.       LOGICAL modified
  125.       CHARACTER naca*(*)
  126. c
  127. c  This routine converts a NACA four didget number into the parameters 
  128. c  describing a foil.
  129. c
  130.       INTEGER ieps,ipmx,itau,iler,im
  131.       INTEGER Ist
  132.       EXTERNAL Ist
  133. c
  134. c     Convert first four digits into integers.
  135. c
  136.       inaca = Ist( 2,naca(3:3),0 )
  137.       ieps  = Ist( 1,naca(1:1),0 )
  138.       ipmx  = Ist( 1,naca(2:2),0 )
  139.       itau  = Ist( 2,naca(3:4),0 )
  140.       ndata(1) = 0.01 * real(ieps)
  141.       ndata(2) =  0.1 * real(ipmx)
  142.       ndata(3) = 0.01 * real(itau)
  143. c
  144. c     Check if a modified section is used.
  145. c
  146.       if (modified) then
  147.      iler = Ist( 1,naca(6:6),6 )
  148.      im   = Ist( 1,naca(7:7),3 )
  149.      ndata(4) = real(iler)
  150.      ndata(5) = 0.1 * real(im)
  151.       endif
  152. c
  153.       return
  154.       END
  155. c
  156. c-----------------------------------------------------------------------------
  157. c
  158.       SUBROUTINE Naca4c( coeff,dcoeff,camber,pmax )
  159. c
  160.       IMPLICIT none
  161.       REAL coeff(3),dcoeff(2),camber,pmax
  162. c
  163. c  This routine calculates the coeff's for the NACA camber equation.
  164. c
  165.       REAL tiny
  166.       PARAMETER( tiny=1.0e-20 )
  167.       REAL inpmax
  168. c
  169. c
  170.       if (pmax.lt.tiny) then
  171.          inpmax = 0.0
  172.       else
  173.          inpmax = 1.0/pmax
  174.       endif
  175. c
  176. c
  177. c
  178.       coeff(1)  = 0.0
  179.       if (camber.lt.tiny) then
  180.      coeff(2) = 0.0
  181.      coeff(3) = 0.0
  182.       else
  183.      coeff(2)  = 2.0*camber*inpmax
  184.      coeff(3)  = (-camber)*(inpmax**2)
  185.       endif
  186.       dcoeff(1) = coeff(2)
  187.       dcoeff(2) = 2.0*coeff(3)
  188. c
  189.       return
  190.       END
  191. c
  192. c-----------------------------------------------------------------------------
  193. c
  194.       SUBROUTINE Naca5( x,npl,npmx,naca,inaca,scle,modified )
  195. c
  196.       IMPLICIT none
  197.       INTEGER npl,npmx,inaca
  198.       REAL x(3,npmx),scle
  199.       CHARACTER naca*(*)
  200.       LOGICAL modified
  201. c
  202. c  This routine calculates the ordinates of a NACA five digit foil.
  203. c
  204.       INTEGER npt,npf
  205.       INTEGER i,j
  206.       REAL ndata(7)
  207.       REAL a(5),d(4),cl,k,m,xi,ct(5)
  208.       REAL dx,dy,dydx,t,y
  209.       REAL Poly
  210.       EXTERNAL Poly
  211. c
  212.       DATA  ct / 0.0 ,-0.126 ,-0.3537 , 0.2843 ,-0.1015 /
  213. c
  214. c
  215.       npt  = npl+1
  216.       npf  = npl/2
  217. c
  218. c     Set the x values according to a cosine distribution.
  219. c
  220.       call cosdn( x,npl,npmx )
  221. c
  222. c     Calculate the parameters describing the foil, from the given NACA
  223. c     number. Then calculate the polynomial coeffs for the camber.
  224. c
  225.       call Naca5b( naca,inaca,ndata,modified )
  226.       if (modified) then
  227.      call Nacam( a,d,ndata(6) )
  228.       endif
  229.       cl = ndata(1)
  230.       k  = ndata(5)
  231.       m  = ndata(4)
  232. c
  233. c     Calculate 'y' values for each value of 'x'.
  234. c
  235.       x(2,1)     = 0.0
  236.       x(2,npt)   = 0.0
  237.       x(2,npf+1) = 0.0
  238.       do i = 2,npf
  239.      j = npt+1-i
  240. c
  241. c        Calculate the camber line.
  242. c
  243.      xi = x(1,i)
  244.      if (xi.le.m) then
  245.         y    =k*cl*(  xi*xi*xi - 3.0*m*xi*xi + m*m*(3.0-m)*xi )/1.8
  246.         dydx =k*cl*( 3.0*xi*xi - 6.0*m*xi    + m*m*(3.0-m)    )/1.8
  247.      else
  248.         y    = k*cl*m*m*m*( 1.0-xi )/1.8
  249.         dydx =-k*cl*m*m*m/1.8
  250.      endif
  251. c
  252. c        Calculate thickness.
  253. c
  254.      if (modified) then
  255.         if (x(1,i).le.ndata(7)) then
  256.            a(1) = a(5)*sqrt( x(1,i) )
  257.            t = ndata(3)*Poly( x(1,i),3,a )
  258.         else
  259.            t = ndata(3)*Poly( (1.0-x(1,i)),3,d )
  260.         endif
  261.      else
  262.         ct(1) = 0.2969 * sqrt( x(1,i) )
  263.         t = 5.0*ndata(3)*Poly( x(1,i),4,ct )
  264.      endif
  265. c
  266. c        Finally, calculate node positions upon the foil surface.
  267. c
  268.      dy = sqrt(( t**2 )/( 1.0+( dydx**2 )))
  269.      dx = dy * dydx
  270. c
  271.      x(1,i) = x(1,i) + dx
  272.      x(1,j) = x(1,i) - 2.0*dx
  273.      x(2,i) = y - dy
  274.      x(2,j) = y + dy
  275.       enddo
  276. c
  277. c     Rescale dimensions if nessisary.
  278. c
  279.       call FoilScale( x,npt,scle )
  280. c
  281.       return
  282.       END
  283. c
  284. c-----------------------------------------------------------------------------
  285. c
  286.       SUBROUTINE Naca5b( naca,inaca,ndata,modified )
  287. c
  288.       IMPLICIT none
  289.       INTEGER inaca
  290.       REAL ndata(7)
  291.       CHARACTER naca*(*)
  292.       LOGICAL modified
  293. c
  294. c  This routine converts a NACA five didgit number into the parameters 
  295. c  describing a foil.
  296. c
  297.       INTEGER icld,ipmt,itau,iler,im
  298.       REAL c1(4),c2(6)
  299.       INTEGER Ist
  300.       REAL Poly
  301.       EXTERNAL Poly,Ist
  302. c
  303.       DATA c1 / 0.0000,   1.0988, 1.3984, 1.8519 /
  304.       DATA c2 / 0.0000,-.0069567, .47172, 16.062,-17.409, 101.24 /
  305. c
  306. c     Convert first four digits into integers.
  307. c
  308.       inaca= Ist( 2,naca(4:5),0 )
  309.       icld = Ist( 1,naca(1:1),0 )
  310.       ipmt = Ist( 2,naca(2:3),0 )
  311.       itau = Ist( 2,naca(4:5),0 )
  312.       ndata(1) =  0.15*real(icld)
  313.       ndata(2) = 0.005*max( real(ipmt),1.0 )
  314.       ndata(3) =  0.01*real(itau)
  315.       ndata(4) = Poly( ndata(2),3,c1 )
  316.       ndata(5) = 1.0/( Poly( ndata(2),5,c2 ) )
  317. c
  318. c     Check if a modified section is used.
  319. c
  320.       if (modified) then
  321.      iler = Ist( 1,naca(7:7),6 )
  322.      im   = Ist( 1,naca(8:8),3 )
  323.      ndata(6) = real(iler)
  324.      ndata(7) = 0.1*real(im)
  325.       endif
  326. c
  327.       return
  328.       END
  329. c
  330. c-----------------------------------------------------------------------------
  331. c
  332.       SUBROUTINE Nacam( a,d,ndata )
  333. c
  334.       IMPLICIT none
  335.       REAL a(0:4),d(0:3),ndata(2)
  336. c
  337. c  This routine calculates the coefficients for the NACA 4/5 digit
  338. c  modified thickness distributions.
  339. c
  340.       REAL mat(3,3),vector(3)
  341.       REAL coeff(4),m,n
  342.       REAL Poly
  343.       EXTERNAL Poly,Crm3
  344. c
  345.       DATA coeff / 1.0090,-0.24048,-2.1786, 15.833 /
  346. c
  347.       m = ndata(2)
  348.       n = 1.0-m
  349. c
  350. c     Calculate coefficients for the aft thickness distribution.
  351. c
  352.       d(0) = 0.0
  353. cc    d(0) = 0.01
  354.       d(1) = Poly( m,3,coeff )
  355.       d(2) = ( 3.0*( 0.5 - d(0) ) - 2.0*d(1)*n )/( n*n )
  356.       d(3) = ( 2.0*( 0.5 - d(0) ) - 1.0*d(1)*n )/(-n*n*n )
  357. c
  358. c     Calculate coefficients for the fore thickness.
  359. c
  360.       a(4) = 0.24742*ndata(1)
  361.       vector(1) = 0.5 - a(4)*sqrt(m)
  362.       vector(2) =-0.5*a(4)/sqrt(m)
  363.       vector(3) = 2.*d(2) + 6.*( 1.0-m )*d(3) + 0.25*a(4)/( m**1.5 )
  364.       mat(1,1) = m
  365.       mat(1,2) = m*m
  366.       mat(1,3) = m*m*m
  367.       mat(2,1) = 1.0
  368.       mat(2,2) = 2.0*m
  369.       mat(2,3) = 3.0*m*m
  370.       mat(3,1) = 0.0
  371.       mat(3,2) = 2.0
  372.       mat(3,3) = 6.0*m
  373.       call Crm3( mat,vector,a(1) )
  374. c
  375.       return
  376.       END
  377.